Tech Arsenal 1

home *** CD-ROM | disk | FTP | other *** search

/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-19 / trjmkr.zip / MATH20.DOC < prev next >

Wrap

Text File | 1992-08-09 | 109KB | 3,607 lines

TRAJECTORY MAKER A Trajectory Simulation Tool ──────────────────────────── ALGORITHM SUPPLEMENT ORBIT HORIZONS 94 Promenade Irvine, California 92715 U.S.A. COPYRIGHT NOTICE This document and software package consisting of the Trajectory Maker and Trajectory Scape computer programs are copyrighted (C) 1991, 1992 by H.B. Reynolds, President, ORBIT HORIZONS. All rights are reserved. No part of this publication may be reproduced, transmitted, transcribed, stored in any retrieval system, or translated into any language by any means without the express written permission of ORBIT HORIZONS, 94 Promenade, Irvine California 92715, USA. FIRST EDITION/FIRST PRINTING July 1992 ii LIMITED WARRANTY THIS ALGORITHM SUPPLEMENT, THE USER'S GUIDE, AND THE PROGRAMS, "TRAJECTORY MAKER" AND "TRAJECTORY SCAPE" ARE SOLD "AS IS", WITHOUT WARRANTY AS TO THEIR PERFORMANCE, OR FITNESS FOR ANY PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE RESULTS AND PERFORMANCE OF THE PROGRAMS IS ASSUMED BY THE USER. HOWEVER, ORBIT HORIZONS WARRANTS THE MAGNETIC DISKETTE(S) ON WHICH THE PROGRAM IS RECORDED TO BE FREE FROM DEFECTS IN MATERIALS AND FAULTY WORKMANSHIP UNDER NORMAL USE FOR A PERIOD OF NINETY DAYS FROM THE DATE OF PURCHASE. IF DURING THIS NINETY DAY PERIOD THE DISKETTE SHOULD BECOME DEFECTIVE, IT MAY BE RETURNED TO ORBIT HORIZONS FOR A REPLACEMENT WITHOUT CHARGE, PROVIDED YOU SEND PROOF OF PURCHASE OF THE PROGRAM. YOUR SOLE AND EXCLUSIVE REMEDY IN THE EVENT OF A DEFECT IS EXPRESSLY LIMITED TO REPLACEMENT OF THE DISKETTE AS PROVIDED ABOVE. IF FAILURE OF A DISKETTE HAS RESULTED FROM ACCIDENT OR ABUSE, ORBIT HORIZONS SHALL HAVE NO RESPONSIBILITY TO REPLACE THE DISKETTE UNDER THE TERMS OF THIS LIMITED WARRANTY. IN NO EVENT SHALL ORBIT HORIZONS BE LIABLE TO YOU FOR ANY DAMAGES, INCLUDING BUT NOT LIMITED TO ANY LOST PROFITS, LOST SAVINGS, OR OTHER CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THESE PROGRAMS AND ACCOMPANYING MATERIALS EVEN IF ORBIT HORIZONS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. SOME STATES DO NOT ALLOW THE LIMITATION OR EXCLUSION OF LIABILITY FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THE ABOVE LIMITATION OR EXCLUSION MAY NOT APPLY TO YOU. THIS WARRANTY GIVES YOU SPECIFIC LEGAL RIGHTS, AND YOU MAY ALSO HAVE OTHER RIGHTS WHICH VARY FROM STATE TO STATE. THIS WARRANTY SHALL BE CONSTRUED, INTERPRETED AND GOVERNED BY LAWS OF THE STATE OF CALIFORNIA. YOU AGREE TO THESE TERMS BY YOUR DECISION TO USE THIS SOFTWARE. iii iii TABLE OF CONTENTS 1 INTRODUCTION ...................................... 1-1 2 EQUATIONS OF MOTION ............................... 2-1 3 GRAVITY PERTURBATION MODEL ........................ 3-1 4 DRAG PERTURBATION MODEL ........................... 4-1 4.1 Density Profile ............................. 4-4 4.2 Wind Profile ............................... 4-14 5 NUMERICAL INTEGRATION SCHEME ...................... 5-1 6 COORDINATE SYSTEMS AND TRANSFORMATIONS ............ 6-1 6.1 Time Considerations ......................... 6-1 6.2 Orbital Elements ............................ 6-2 6.3 Geospherical Inertial Coordinates ........... 6-6 6.4 Earth-Fixed-Geocentric Coordinates .......... 6-6 6.5 Local Topocentric Coordinates ............... 6-7 6.6 Common Radar Coordinates ................... 6-10 6.7 Geodetic Coordinates ....................... 6-10 7 BIBLIOGRAPHY ....................................... 7-1 LIST OF TABLES Table Page 3-1 WGS-72 Earth Constants ........................ 3-1 3-2 Smithsonian Zonal Harmonic Coefficients ....... 3-2 4-1 Drag Coefficient vs. Mach Number .............. 4-3 4-2 Molecular Scale Temperature ................... 4-4 4-3 Atmospheric Density Profile ................... 4-6 5-1 Shanks 8-12 Numerical Integration Constants ... 5-3 iv LIST OF FIGURES Figure Page 3-1 Magnitude of Vector Differences Between the 23rd and 2nd, 3rd, 4th and 5th Order Gravity Models ......................................... 3-7 3-2 Magnitude of Vector Differences Between the 23rd and 6th, 7th, 8th and 9th Order Gravity Models ......................................... 3-8 3-3 Magnitude of Vector Differences Between the 23rd and 10th, 11th, 12th and 13th Order Gravity Models ................................. 3-9 3-4 Magnitude of Vector Differences Between the 23rd and 14th, 15th, 16th and 17th Order Gravity Models ........................................ 3-10 3-5 Magnitude of Vector Differences Between the 23rd and 18th, 19th, 20th and 21st Order Gravity Models ................................ 3-11 3-6 Magnitude of Vector Differences Between the 23rd and 19th, 20th, 21st and 22nd Order Gravity Models ................................ 3-12 3-7 Magnitude of Vector Differences Between the 23rd and 2nd, 4th, 6th and 8th Order (Even) Gravity Models ................................ 3-13 3-8 Magnitude of Vector Differences Between the 23rd and 3rd, 5th, 7th and 9th Order (Odd) Gravity Models ................................ 3-14 3-9 Magnitude of Vector Differences Between the 23rd and 2nd, 6th, 10th and 14th Order Gravity Models ................................ 3-15 3-10 Magnitude of Vector Differences Between the 23rd and 5th, 10th, 15th and 20th Order Gravity Models ................................ 3-16 4-1 1976 U. S. Standard Atmosphere Density Profile ........................................ 4-7 4-2 Abbreviated 1976 U. S. Standard Atmosphere Density Profile ................................ 4-8 v LIST OF FIGURES Figure Page 4-3 Trajectory Maker Density Profile Model ......... 4-9 4-4 Cubic Spline Interpolation Error for 21 and 48 Point Tables ............................... 4-10 4-5 Cubic Spline Interpolation Error for 21 and 48 Point Tables ............................... 4-11 4-6 Cubic Spline and Linear Interpolation Error for 48 Point Table ............................ 4-12 4-7 Cubic Spline Interpolation Error for 48 Point Table ......................................... 4-13 4-8 East/North Wind velocity Components (Example 1) ................................... 4-15 4-9 East/North Wind velocity Components (Example 2) ................................... 4-16 5-1 Numerical Integration Error .................... 5-6 vi Trajectory Maker Introduction Algorithm Supplement 1-1 ────────────────────────────────────────────────────────── 1. INTRODUCTION The Trajectory Maker/Trajectory Scape software system is a state of the art simulation and design tool for use on IBM compatible personal computers. Within minutes it predicts precision trajectory ephemerides and displays their ground traces on world map projections. The shapes of various orbit ground traces can be quite bizarre, unexpected and difficult to visualize because of two motion processes; namely the Earth's rotation about its axis, coupled with the motion of the orbiting object. With Trajectory Maker, visualization of this complex interaction is a breeze. But Trajectory Maker gives much more than this graphics' visualization capability. It also creates high-precision numerical output data, in several formats, for those wanting more comprehensive analysis capabilities. Major features of the software and how to use it are described in detail in the User's Guide. The User's Guide also contains some historical background material, a glossary of terms, definitions and illustrations to aid the user with the external interface to the system. Transparent to the user is an underlying foundation of sophisticated mathematical algorithms and techniques that are employed by the system in order to perform the relevant numerical computations, in an accurate and expedient manner. The user of the system need not necessarily be concerned with its inner workings; however, an understanding of theory is sometimes helpful in understanding its practice. This paper presents a detailed description of mathematical algorithms and techniques used by the Trajectory Maker/ Trajectory Scape software system. In addition to presenting the mathematical core of the software system, the merits and limitations of the algo- rithms are discussed. Illustrations are provided that will aid in selecting various aspects of the models, such as truncation error, integration step-size, etc., for the particular application under consideration. The equations of motion that govern the trajectories of Trajectory Maker Introduction Algorithm Supplement 1-2 ────────────────────────────────────────────────────────── Earth orbiting satellites and ballistic missiles are summarized. The force models include perturbation effects due to Earth oblateness, aerodynamic drag and atmospheric winds. These forces are developed and investigated. Explicit closed-form expressions for the gravitational acceleration components are derived for an axi-symmetric oblate Earth model. The model contains an arbitrary number of zonal harmonic coefficients. Effects of various orders of the coefficients on certain mission ephemerides are investigated. The aerodynamic drag force is defined, the atmospheric density profile is modeled as a function of altitude and its relative error is analyzed. An innovative algorithm for computing altitude is presented that is very fast, yet produces highly accurate values. A numerical scheme for integrating the equations of motion is presented, relevant coordinate systems are defined and algorithms for transforming the various input and output parameters are summarized. All physical constants used by the software are included. The algorithms presented herein support the trajectory mechanics functions in the program source code with examples that demonstrate their integrity. A bibliography is provided for algorithm extensions, map projection related processes and routine mathematical processes. Trajectory Maker Equations of Motion Algorithm Supplement 2-1 ────────────────────────────────────────────────────────── 2. EQUATIONS OF MOTION The equations of motion for propagating trajectories are modeled in rectangular, geocentric inertial coordinates with fundamental plane coplanar with the Earth's equator and principle axis in the direction of the vernal equinox. Let x, y, z denote these coordinates with z-axis coinci- dent with the Earth's spin-axis, positive northwards; the x-axis in the equatorial plane, positive in the direction of the vernal equinox; and the y-axis lying in the equatorial plane in such a manner that the system is right-handed. Further, let i, j, k be unit vectors along the x, y, z axes, respectively. Then, the position of a point mass m may be represented in vector form r = xi + yj + zk and according to Newton's second law, the differential equations of motion may be written d2r m ─── = Σ F dt2 where Σ F is the sum of the external forces acting on the mass m. In this paper the forces due to the Earth's gravitational potential function and its atmosphere are treated. Letting G and D denote these forces per unit mass, we may write d2r ─── = G + D dt2 By substitution of variables, this equation may be decomposed into the following first order differential equations, dr ── = v dt Trajectory Maker Equations of Motion Algorithm Supplement 2-2 ────────────────────────────────────────────────────────── dv ── = G + D dt where, using Newtonian differentiation, the velocity vector is . . . v = xi + yj + zk Expressed in component form, these two vector differential equations become dx . dy . dz . ── = x , ── = y , ── = z dt dt dt . . . dx dy dz ── = Gx + Dx , ── = Gy + Dy , ── = Gz + Dz dt dt dt A solution to these six differential equations may be obtained, once the components of the gravity and drag acceleration components are expressed as known functions of x, y, z. A time series solution is called an ephem- eris. The actual path of the mass m is its trajectory. Trajectory Maker Gravity Model Algorithm Supplement 3-1 ────────────────────────────────────────────────────────── 3. GRAVITY PERTURBATION MODEL We suppose that the Earth's figure is an oblate spheroid, and let P be point at a distance r from its center of mass and let φ' be the declination of P. Then, the potential at P may be expressed as a series of spherical harmonics of the form µ ┌─ ∞ ─┐ V = - ── │ 1 - Σ Jn ( R/r)n Pn(sin φ') │ r └─ n=2 ─┘ where µ is the Earth's gravitational mass constant, R its equatorial radius, Jn the zonal harmonic coefficients and Pn the associated Legendre polynomials of order n. This form of the Earth's potential function is known as the Vinti potential function [1]. The gravitational mass constant and Earth radii used in this paper are from the WGS-72 constants [3] shown in Table 3-1. Table 3-1. WGS-72 Earth Constants ╒═══╤═══════════════════════╤═════════════════════╕ │ │ English │ Metric │ ├───┼───────────────────────┼─────────────────────┤ │ µ │ 62750.16789 nm3/sec2 │ 398600.50 km3/sec2 │ ├───┼───────────────────────┼─────────────────────┤ │ R │ 3443.91739 nm │ 6378.135 km │ ├───┼───────────────────────┴─────────────────────┤ │ Ω │ 0.7292115147 (10)-4 rad/sec │ ├───┼─────────────────────────────────────────────┤ │ f │ 1/298.26 Dimensionless │ └───┴─────────────────────────────────────────────┘ The harmonic coefficients are the Smithsonian 1973 II coefficients [2] tabulated in Table 3-2. Trajectory Maker Gravity Model Algorithm Supplement 3-2 ────────────────────────────────────────────────────────── Table 3-2. Smithsonian Zonal Harmonic Coefficients ╒════════════════════════╦════════════════════════╕ │ J2 = 1082.636e-6 ║ J13 = -0.336e-6 │ │ J3 = -2.540e-6 ║ J14 = 0.101e-6 │ │ J4 = -1.619e-6 ║ J15 = 0.104e-6 │ │ J5 = -0.230e-6 ║ J16 = 0.043e-6 │ │ J6 = 0.552e-6 ║ J17 = -0.227e-6 │ │ J7 = -0.345e-6 ║ J18 = -0.077e-6 │ │ J8 = -0.204e-6 ║ J19 = 0.083e-6 │ │ J9 = -0.162e-6 ║ J20 = -0.108e-6 │ │ J10 = -0.232e-6 ║ J21 = -0.070e-6 │ │ J11 = 0.317e-6 ║ J22 = 0.075e-6 │ │ J12 = -0.196e-6 ║ J23 = 0.111e-6 │ └────────────────────────╨────────────────────────┘ Now, the distance from the Earth's center to the mass m, coincident with the point P, is the magnitude of r, given by r = √(x2 + y2 + z2) and the sine of its declination is sin φ' = z/r Since the gravitational force is conservative, the net gravitational acceleration G acting on the mass is the negative gradient of V G = - V which can be decomposed, by differentiating, into the form G = -(µ/r3)r - U The first term in this expression represents the dominant attraction the Earth would exhibit if it were a sphere comprised of homogeneous concentric shells. This term accounts for the non-perturbed conic motion of m. The function U is the disturbing function whose gradient represents the gravity perturbation acceleration which causes deviations from the conic motion of m. This function is given by Trajectory Maker Gravity Model Algorithm Supplement 3-3 ────────────────────────────────────────────────────────── ∞ U = µ Σ Jn Rn r-(n+1) Pn(sin φ') n=2 It follows from the gradient operator that the components of the gravitational perturbation acceleration are δGx ≡ - U/ x = (µx)/r3 Σ Jn(R/r)n fn(φ') δGy ≡ - U/ y = (µy)/r3 Σ Jn(R/r)n fn(φ') δGz ≡ - U/ z = (µz)/r3 Σ Jn(R/r)n gn(φ') where fn(sin φ') = (n + 1)Pn(sin φ') + sin φ Pn'(sin φ') gn(sin φ') = (n + 1)Pn(sin φ') - cos φ cot φ Pn'(sin φ') Whittaker and Watson [5] derive the following four recurrence formulas for the complex variable u: Pn+1'(u) - uPn'(u) = (n + 1)Pn(u) (1) (u2 - 1)Pn'(u) = nuPn(u) - nPn-1(u) (2) (n + 1)Pn+1(u) - (2n + 1)uPn(u) + nPn-1(u) = 0 (3) Pn+1'(u) - Pn-1'(u) = (2n + 1)Pn(u) (4) Using the first recurrence formula yields fn(sin φ') = Pn+1'(sin φ') and the second and third formulas yield gn(sin φ') = (n + 1)Pn+1(sin φ')/sin φ' Therefore, the components of the perturbation acceleration become δGx = (µx)/r3 Σ Jn(R/r)n Pn+1'(sin φ') δGy = (µy)/r3 Σ Jn(R/r)n Pn+1'(sin φ') δGz = (µ)/r2 Σ Jn(R/r)n(n + 1)Pn+1(sin φ') In order to numerically evaluate these expressions it is, Trajectory Maker Gravity Model Algorithm Supplement 3-4 ────────────────────────────────────────────────────────── of course, necessary to limit the series' to a finite number of terms. Suppose, then, that the degree of the polynomials in R/r is specified, say M ≥ 2. Then, with a shift of indices we have M-1 δGx = (µx)/r3 (R/r)2 Σ ai(R/r)i-1 i=1 M-1 δGy = (µy)/r3 (R/r)2 Σ ai(R/r)i-1 i=1 M-1 δGz = (µ)/r2 (R/r)2 Σ bi(R/r)i-1 i=1 where ai = Ji+1 Pi+2'(sin φ') bi = Ji+1(i + 2)Pi+2(sin φ') The Legendre polynomials and their derivatives for n ≥ 3 may be respectively computed from the third and fourth recurrence formulas once it is noted that [5] P1(u) = u P2(u) = 1/2 (3u2 - 1) and consequently P1'(u) = 1 P2'(u) = 3u In order to illustrate the effects of various combinations of gravity harmonic coefficients, we consider a nominal Defense Meteorological Satellite Program mission. Its state vector, at orbit insertion, is represented in the following Earth-centered-inertial coordinates: Trajectory Maker Gravity Model Algorithm Supplement 3-5 ────────────────────────────────────────────────────────── . x = 442.151588 nm x = 0.512255 nm/sec . y = 1387.396376 nm y = -3.732104 nm/sec . z = -3611.173591 nm z = -1.373147 nm/sec The satellite's near circular orbit is sun-synchronous with altitude of about 450 nm. The orbit ephemerides were propagated from orbit insertion for a duration of 25000 seconds, in increments of 100 sec. This step-size is optimal, in this case, since any smaller step does not significantly improve the solution. The trajectory ephemeris spans eight orbits. Various gravity models with reduced orders (less coeffi- cients) were compared with the model of highest order having the 22 coefficients listed in Table 2 above. The state vector of the highest order model was compared with the lower order models by differencing their respec- tive position vector components, and computing the magnitude of the resulting difference vector. In the figures the difference between a trajectory generated with a gravity model of order M and one of order N is indicated in the legend by the symbol DJMJN. The first group of the figures, namely, Figures 3-1, 3-2, 3-3, 3-4, 3-5, 3-6, illustrate differences in ephemerides, sequentially beginning with the J2 coefficient and ending with coefficients up to J23. It is seen that this family of curves has a secular growth over the time span, while exhibiting a tendency towards convergence; however, that is not to say that any given order gravity model produces more accurate trajectories than a lower order model, as for example, comparison of the sixth and seventh order models in Figure 3-2. Figure 3-3 illustrates the beginning of the loss of data integrity as round-off accumulates, becoming more pronounced in Figures 3-4 and 3-5. The scales in the above group of figures vary across the figures by two orders of magnitude. In order to illus- trate variations across scales, the next group of figures, namely, Figures 3-7, 3-8, 3-9, and 3-10 were composed. Not all differences are illustrated, but enough to get the idea. Trajectory Maker Gravity Model Algorithm Supplement 3-6 ────────────────────────────────────────────────────────── The gravitational effects on other orbit parameters is given in detail in Reynolds [4]. Trajectory Maker Gravity Model Algorithm Supplement 3-7 ────────────────────────────────────────────────────────── Figure 3-1. Magnitude of Vector Differences Between the 23rd and 2nd, 3rd, 4th and 5th Order Gravity Models. Trajectory Maker Gravity Model Algorithm Supplement 3-8 ────────────────────────────────────────────────────────── Figure 3-2. Magnitude of Vector Differences Between the 23rd and 6th, 7th, 8th and 9th Order Gravity Models. Trajectory Maker Gravity Model Algorithm Supplement 3-9 ────────────────────────────────────────────────────────── Figure 3-3. Magnitude of Vector Differences Between the 23rd and 10th, 11th, 12th and 13th Order Gravity Models. Trajectory Maker Gravity Model Algorithm Supplement 3-10 ────────────────────────────────────────────────────────── Figure 3-4. Magnitude of Vector Differences Between the 23rd and 14th, 15th, 16th and 17th Order Gravity Models. Trajectory Maker Gravity Model Algorithm Supplement 3-11 ────────────────────────────────────────────────────────── Figure 3-5. Magnitude of Vector Differences Between the 23rd and 18th, 19th, 20th and 21st Order Gravity Models. Trajectory Maker Gravity Model Algorithm Supplement 3-12 ────────────────────────────────────────────────────────── Figure 3-6. Magnitude of Vector Differences Between the 23rd and 19th, 20th, 21st and 22nd Order Gravity Models. Trajectory Maker Gravity Model Algorithm Supplement 3-13 ────────────────────────────────────────────────────────── Figure 3-7. Magnitude of Vector Differences Between the 23rd and 2nd, 4th, 6th and 8th Order (Even) Gravity Models. Trajectory Maker Gravity Model Algorithm Supplement 3-14 ────────────────────────────────────────────────────────── Figure 3-8. Magnitude of Vector Differences Between the 23rd and 3rd, 5th, 7th and 9th Order (Odd) Gravity Models. Trajectory Maker Gravity Model Algorithm Supplement 3-15 ────────────────────────────────────────────────────────── Figure 3-9. Magnitude of Vector Differences Between the 23rd and 2nd, 6th, 10th and 14th Order Gravity Models. Trajectory Maker Gravity Model Algorithm Supplement 3-16 ────────────────────────────────────────────────────────── Figure 3-10. Magnitude of Vector Differences Between the 23rd and 5th, 10th, 15th and 20th Order Gravity Models. Trajectory Maker Gravity Model Algorithm Supplement 3-17 ────────────────────────────────────────────────────────── References 1. Anon., Space Flight Handbooks, Volume I, NASA SP-33 Part 1, Office of Scientific and Technical Inform- ation, NASA, Washington D. C., 1963. 2. Gaposchkin, F. M., ed., 1973 Smithsonian Standard Earth (III), Smithsonian Astrophysical Observatory, Special Report 353, Cambridge, Ma., November 28, 1973. 3. Hieb, H., Western Test Range Geodetic Coordinate Manual, Part 1, Federal Electric Corporation, Feb. 1980. 4. Reynolds, H. B., High-Order Gravitational Perturbation Effects On DMSP Extended Mission Ephemerides, OAO Corporation, Technical Operating Report, El Segundo, Ca., Jan. 3, 1985. 5. Whittaker, F. and G. Watson, A Course in Modern Analysis, Cambridge University Press, Fourth Ed. 1963. Trajectory Maker Drag Model Algorithm Supplement 4-1 ────────────────────────────────────────────────────────── 4. DRAG PERTURBATION MODEL The aerodynamic drag force acting on an object varies jointly with the atmospheric density σ, the object drag coefficient CD, the projected cross sectional area A of the object (normal to the velocity vector), the square of its speed relative to the rotating atmosphere va, and inversely with its mass m. The drag force acts in a direction opposite that of the relative velocity vector. In vector form this force per unit mass may be written CD A D = - σ ──── va va 2m justifiably, since the object imparts a velocity of order va to a mass of air Aσva in each second; the rate of change of momentum, equal to the drag force on the object, is thus Aσva². The empirical drag coefficient is based on such things as the object's shape, altitude and Mach number. Taking into account local wind conditions with the vector q, the velocity relative to the rotating air mass is va = v - Ω x r - q where Ω = Ωk is the Earth's rotation rate vector having magnitude 0.7292115147 (10)-4 rad/sec. It is usual to define a ballistic coefficient ß by m ß = ──── CD A and write the drag force as 1 D = - ─── σ va va 2 ß which in component form, after expanding the above vector product, becomes Trajectory Maker Drag Model Algorithm Supplement 4-2 ────────────────────────────────────────────────────────── . Dx = -σva[(x + Ωy) - q1]/(2ß) . Dy = -σva[(y - Ωx) - q2]/(2ß) . Dz = -σva[ z - q3 ]/(2ß) where . . . va = √{[(x + Ωy) - q1]2 + [(y - Ωx) - q2]2 + (z - q3)2} and q1, q2, q3 denote the wind velocity components relative to the rotating atmosphere. Ballistic coefficients have an intuitive appeal since their units, measured in kg/m², resemble density units. When working with English units, the ballistic coefficient is usually expressed in mean sea level weight rather than mass, that is lb/ft². The following are typical examples that may be encountered: o A space based radar antenna: 0.144, to 0.448 lb/ft² o Space Shuttle Discovery: 42 lb/ft² o Ballistic missiles: 25, 75, 2150 lb/ft² When computing a ballistic coefficients for orbiting objects, the parameter with the most uncertainty is usually the drag coefficient. Although this is not always the case. For example, Reynolds [2] considers a space based radar system having a rotating antenna whose shape is an arbitrary surface of revolution. Wertz [5] suggests that in the absence of aposteriori knowledge of the drag coefficient that a value of 2.0 be used. Waison [4] reports that it should be consistent with the value used to derive the atmospheric model, and that for the 1959 ARDC model one should use 2.0; for the 1964 Jacchia model, a value of 2.2 should be used. For reentry vehicles, the drag coefficient may be provided as a function of Mach number M, as in the Table 4-1, together with the vehicle's mass and effective area. Trajectory Maker Drag Model Algorithm Supplement 4-3 ────────────────────────────────────────────────────────── Table 4-1. Drag Coefficient vs. Mach Number Mass = 35.8365 slugs, Area = 19.8 ft² ┌────╥─────┬─────┬─────┬─────┬─────┬──────┬─────┬─────┐ │ M ║ 0 │ .25 │ .50 │ .75 │ 1.0 │ 1.25 │ 1.5 │ 100 │ ├────╫─────┼─────┼─────┼─────┼─────┼──────┼─────┼─────┤ │ CD ║ .38 │ .40 │ .44 │ .55 │ .72 │ .76 │ .77 │ .77 │ └────╨─────┴─────┴─────┴─────┴─────┴──────┴─────┴─────┘ Mach number is function of the speed of sound, which in turn is a function of temperature. The expression adopted for the speed of sound in Reference [1] is Cs = √( R* TM/M0) where is the ratio of specific heat of air at constant pressure to that at constant volume, and taken to be 1.4 (dimensionless) exactly, R* the gas constant taken to be 8314.32 N.m/(kmol.K), M0 the sea-level mean molecular weight of the mixture of gasses taken to be 28.9644 kg/kmol, and TM is the molecular scale temperature, modeled as a linear, segmental arc. Putting all these things together results in Cs = 20.04680276 √TM (meters/sec) and the Mach number is M = va/Cs The concept of speed of sound progressively looses its range of applicability at high altitudes and for this reason the values of speed of sound are not computed for heights above 86 km. Table 4-2 is extracted from Table 9 of the U. S. Standard 1976 abbreviated tables. Trajectory Maker Drag Model Algorithm Supplement 4-4 ────────────────────────────────────────────────────────── Table 4-2 Molecular Scale Temperature ╒═══════════════╤══════════════════╕ │ Geometric │ Molecular Scale │ │ Altitude (km) │ Temperature (K) │ ├───────────────┼──────────────────┤ │ 0.0000 │ 288.150 │ │ 11.0190 │ 216.650 │ │ 20.0631 │ 216.650 │ │ 32.1619 │ 288.650 │ │ 47.3500 │ 270.650 │ │ 51.4124 │ 270.650 │ │ 71.8019 │ 214.650 │ │ 86.0000 │ 186.946 │ └───────────────┴──────────────────┘ Since the U. S. STandard models the temperature as linear segments, it is only necessary to interpolate linearly in this table to obtain TM for arbitrary altitude, and substitute the result into the above equation to compute the speed of sound. 4.1 Atmospheric Density Profile Another complication is that density varies in a fairly irregular manner with factors such as solar activity, yearly seasons, time of day, geographical location and altitude. Fortunately, extensive research has been performed by various agencies to define models for the principle atmospheric parameters. The 1976 U. S. Standard Atmosphere [1] is an idealized, steady state representation of the Earth's atmosphere from its surface to 1000 km, as it is assumed to exist in a period of moderate solar activity. It consists of compre- hensive tables (approximately 1000 points) of various atmospheric parameters, including density, tabulated as a function of altitude. The density profile is illustrated in Figure 4-1 on a common logarithm scale. Since 1000 pairs is a lot of data to handle, subsets were extracted and modeled with the goal of obtaining a rela- tively small subset that would yield sufficient accuracy. Trajectory Maker Drag Model Algorithm Supplement 4-5 ────────────────────────────────────────────────────────── In effort to obtain the smoothest possible representation for density, cubic spline functions were used to model the common logarithm of the density profile, with the realiza- tion that, during orbit propagation, density would necessarily be obtained from the logarithm profile by exponentiation. Several subsets were extracted from the Standard Atmo- sphere, starting with the abbreviated tables, consisting of 21 pairs, from the Standard Atmosphere [1]. These points are illustrated in Figure 4-2. The final set selected for the model, herein, was obtained by judiciously selecting points from the idealized model until the percent relative error in density was reduced to less than 1%. The resulting set consisted of 48 altitude- density pairs illustrated in Figure 4-3. Note that these data are not equally spaced. However, besides being very fast, Thompson's algorithm [3] accounts for this situa- tion. Figure 4-4 illustrates the percent relative error obtained when using these two density profiles to model the idealized atmosphere. The dashed curve shows that in the critical altitude region (less than 200 km) for reentry trajectories, the 21 point set produces errors as high as 6%. On the other hand, for the 48 point set the errors are less than 1%. This is shown more clearly in Figure 4-5. Figure 4-6 compares linear interpolation (dashed curve) with spline interpolation (solid curve). It is seen that linear interpolation produces errors as high as 4%. It should be emphasized that these error curves are density errors rather than the logarithm of the density errors. Figure 4-7 isolates the error for the 48-point set on a magnified scale. The figure shows that the error to be less than 1% on the entire altitude range. It is this set that comprises the density profile used herein. The actual values for this set are tabulated in Table 4-3. This study also found that the endpoint derivatives played a significant role towards achieving the 1% accuracy criterion. Thus, rather than computing the derivatives with values from Table 4-3, the following values were computed using a finer resolution from data in the comprehensive table: Trajectory Maker Drag Model Algorithm Supplement 4-6 ────────────────────────────────────────────────────────── (log σ(0))' = -0.041934 (log σ(1000))' = -0.001834 where we recall that the splines were fit to the common logarithm of the density values. Table 4-3. Atmospheric Density Profile ╒════════════════════════╦═════════════════════════╕ │ Altitude Density ║ Altitude Density │ │ (km) (kg/m3) ║ (km) (kg/m3) │ ├────────────────────────╫─────────────────────────┤ │ 0 1.2250 ║ 90 3.416e-6 │ │ 2 1.0066 ║ 100 5.604e-7 │ │ 4 8.1935e-1 ║ 110 9.708e-8 │ │ 6 6.6011e-1 ║ 120 2.222e-8 │ │ 8 5.2579e-1 ║ 130 8.152e-9 │ │ 10 4.1351e-1 ║ 140 3.831e-9 │ │ 12 3.1194e-1 ║ 150 2.076e-9 │ │ 14 2.2786e-1 ║ 160 1.233e-9 │ │ 16 1.6647e-1 ║ 170 7.815e-10 │ │ 18 1.2165e-1 ║ 180 5.194e-10 │ │ 20 8.8910e-2 ║ 190 3.581e-10 │ │ 25 4.0084e-2 ║ 200 2.541e-10 │ │ 30 1.8410e-2 ║ 220 1.367e-10 │ │ 35 8.4634e-3 ║ 240 7.858e-11 │ │ 40 3.9957e-3 ║ 260 4.742e-11 │ │ 45 1.9663e-3 ║ 280 2.971e-11 │ │ 50 1.0269e-3 ║ 300 1.916e-11 │ │ 45 1.9663e-3 ║ 400 2.802e-12 │ │ 60 3.0968e-4 ║ 500 5.215e-13 │ │ 65 1.6321e-4 ║ 600 1.137e-13 │ │ 70 8.2829e-5 ║ 700 3.069e-14 │ │ 75 3.9921e-5 ║ 800 1.136e-14 │ │ 80 1.8458e-5 ║ 900 5.759e-15 │ │ 85 8.2196e-6 ║ 1000 3.561e-15 │ └────────────────────────╨─────────────────────────┘ Trajectory Maker Drag Model Algorithm Supplement 4-7 ────────────────────────────────────────────────────────── Fig. 4-1. 1976 U. S. Standard Atmosphere Density Profile. Trajectory Maker Drag Model Algorithm Supplement 4-8 ────────────────────────────────────────────────────────── Fig. 4-2. Abbreviated 1976 U. S. Standard Atmosphere Density Profile. Trajectory Maker Drag Model Algorithm Supplement 4-9 ────────────────────────────────────────────────────────── Fig. 4-3. Trajectory Maker Density Profile Model. Trajectory Maker Drag Model Algorithm Supplement 4-10 ────────────────────────────────────────────────────────── Fig. 4-4. Cubic Spline Interpolation Error for the 21 and 48 Point Tables. Trajectory Maker Drag Model Algorithm Supplement 4-11 ────────────────────────────────────────────────────────── Fig. 4-5. Cubic Spline Interpolation Error for the 21 and 48 Point Tables. Trajectory Maker Drag Model Algorithm Supplement 4-12 ────────────────────────────────────────────────────────── Fig. 4-6. Cubic Spline and Linear Interpolation Error for the 48 Point Table. Trajectory Maker Drag Model Algorithm Supplement 4-13 ────────────────────────────────────────────────────────── Fig. 4-7. Cubic Spline Interpolation Error for the 48 Point Table. Trajectory Maker Drag Model Algorithm Supplement 4-14 ────────────────────────────────────────────────────────── 4.2 Wind Profile Local atmospheric winds are typically supplied in the form of a table which contains the wind velocity magnitude q(h) and direction Θ(h) as functions of altitude. The direc- tion of the wind is defined as the true bearing from which it is blowing. The east, north and up components of the wind velocity are, respectively, qe = -q(h) sin Θ(h) qn = -q(h) cos Θ(h) qu = 0 For the purpose of illustration, a couple of typical wind profiles from the Lompoc region are illustrated in Figures 4-8 and 4-9. The symbols represent observed wind velocities that were mapped into east and north components by the above equations. Cubic spline functions were used to model these data and several values were computed and plotted at intermediate values of altitude. Two things are noteworthy - that splines provide a smooth represent- ation of the data, and that Lompoc is a windy place. The above components can be transformed into Earth-fixed- geocentric coordinates with the aid of a rotation matrix R which shall be defined, subsequently, when coordinate systems and transformations are treated. At that time, an innovative algorithm for computing altitude will also be developed. For now, it suffices to write ┌─ ─┐ ┌─ ─┐ │ q1 │ │ qe │ │ q2 │ = R │ qn │ │ q3 │ │ qu │ └─ ─┘ └─ ─┘ Trajectory Maker Drag Model Algorithm Supplement 4-15 ────────────────────────────────────────────────────────── Fig. 4-8. East-North Wind Velocity Components, Example 1. Trajectory Maker Drag Model Algorithm Supplement 4-16 ────────────────────────────────────────────────────────── Fig. 4-9. East-North Wind Velocity Components, Example 2. Trajectory Maker Drag Model Algorithm Supplement 4-17 ────────────────────────────────────────────────────────── References 1. Menzner, et. al., Defining Constants, Equations, and Abbreviated Tables of the 1975 U. S. Standard Atmosphere, National Aeronautics and Space Administration, Technical Report NASA TR R-459, Washington D. C., May 1976. 2. Reynolds, H. B., Atmospheric Drag Effects on a Space Based Radar Antenna, TRW Technical Memorandum, January, 1991. 3. Thompson, R. F., Spline Interpolation on a Digital Computer, Goddard Space Flight Center, Technical Report X-692-70-261, Greenbelt, Md., July 1970. 4. Waison, G. R., Prediction of Orbital Decay Rates Due to Atmospheric Drag, TRW Program Technical Report, 3150-6020-RO-000, January 7, 1966. 5. Wertz, J. R., ed., Spacecraft Attitude Determination and Control, D. Reidel Publishing Company, Boston, 1985. Trajectory Maker Integration Scheme Algorithm Supplement 5-1 ────────────────────────────────────────────────────────── 5. NUMERICAL INTEGRATION SCHEME Substituting the atmospheric drag and gravitational perturbation acceleration components, into the component form of the first order vector differential equations of motion, results in the following explicit system of six first-order ordinary differential equations: dx . ── = x dt dy . ── = y dt dz . ── = z dt . dx µx ── = - ── + δGx + Dx dt r3 . dy µy ── = - ── + δGy + Dy dt r3 . dz µz ── = - ── + δGz + Dz dt r3 If the perturbation components are set to zero, the resulting equations are clearly recognizable as describing classical two-body motion. By defining the system state as the vector . . . x = [ x y z x y z ]T where the superscript T denotes matrix transposition, the differential equations may now be brought into a computational convenient vector form dx ── = f(x) dt Trajectory Maker Integration Scheme Algorithm Supplement 5-2 ────────────────────────────────────────────────────────── where f represents the vector-valued function of the state vector defined by the right-hand-side of the differential equations. The differential equations of motion comprise an initial value system, which may be written in more generally in the vector form including time, dx ── = f(t, x), x(t0) = x0 dt . . . where x = (x1, x2, x3, x4, x5, x6) ≡ (x, y, z, x, y, z). A numerical scheme for integrating initial value systems that is of very high-precision, generally stable and relatively easy to implement is Shanks' explicit 8-12 formulation for solutions to differential equations by evaluations of functions [1]. This method is based on the classical Runge-Kutta formulas found in most textbooks on numerical analysis. The 8-12 means that the order of the integration scheme is eight (that is, the truncation error is order nine) and twelve evaluations (that is, stages) of the right-hand side of the differential equations are performed. This scheme is a fixed step integrator which is preferred for trajectory animation graphics. For example, an orbiting object may move farther in its present time step than its preceding time step. To an observer, the object appears to move faster. An advantage of employing a high-order scheme is that accuracy may be maintained in considerable fewer steps than the classical formulations. Given the state x = x(t) of the system at time t, then the state of the system at t + δt is given by x(t + δt) = x(t) + δx where δx ≈ ( 41 k1 + 216 k6 + 272 k7 + 27 k8 + 27 k9 + 36 k10 + 180 k11 + 41 k12 ) / 840 + O(δt9) Trajectory Maker Integration Scheme Algorithm Supplement 5-3 ────────────────────────────────────────────────────────── Using the Einstein summation convention with σ = 1, ..., 12 and µ ≤ σ, the coefficients may be defined by kσ = f(t + tσ δt, x + sσµ kµ ) δt The non-zero constants tσ and sσµ appearing in this expression are defined in Table 5-1 below. Table 5-1. Shanks' 8-12 Numerical Integration Constants. ╒════════════════════════════════════════════════════════╕ │ t2 = 1/9 t3 = 1/6 t4 = 1/4 t5 = 1/10 │ │ t6 = 1/6 t7 = 1/2 t8 = 2/3 t9 = 1/3 │ │ t10= 5/6 t11= 5/6 t12= 1 │ ├────────────────────────────────────────────────────────┤ │ s11 = 1/9 │ │ s21 = 1/24 s22 = 1/8 │ │ s31 = 1/16 s33 = 3/16 │ │ s41 = 29/500 s43 = 33/500 s44 = -3/125 │ │ s51 = 11/324 s54 = 1/243 s55 = 125/972 │ ├────────────────────────────────────────────────────────┤ │ s61 = -7/12 s64 = 19/9 s65 = 125/36 │ │ s66 = -9/2 │ ├────────────────────────────────────────────────────────┤ │ s71 = -10/81 s74 = -32/243 s75 = 125/243 │ │ s77 = 11/27 │ ├────────────────────────────────────────────────────────┤ │ s81 = 1175/324 s84 = -32/3 s85 = -3125/162 │ │ s86 = 26 s87 = 121/162 s88 = -1/12 │ ├────────────────────────────────────────────────────────┤ │ s91 = 293/324 s94 = -71/27 s95 = -1375/324 │ │ s96 = 51/9 s97 = -59/162 s98 = 1/2 │ │ s99 = 1 │ ├────────────────────────────────────────────────────────┤ │ s101 = 1303/1620 s104 = -71/27 s105 = -1375/324 │ │ s106 = 37/6 s107 = 103/162 s1010= 1/10 │ ├────────────────────────────────────────────────────────┤ │ s111 = -955/492 s114 = 2560/369 s115 = 8125/738 │ │ s116 = -612/41 s117 = 7/82 s118 = -27/164 │ │ s119 = -18/41 s1110= -12/41 s1111= 30/41 │ └────────────────────────────────────────────────────────┘ Then, given the initial state x0 of the system at time t0, the time increment δt and the final time tf, the solution may be generated (or propagated) with the following recursion scheme: Trajectory Maker Integration Scheme Algorithm Supplement 5-4 ────────────────────────────────────────────────────────── x(tm) = x(t0 + m δt); m = 0, 1, ..., (tf - t0)/δt A numerical scheme that integrates the differential equations of motion directly, such as this one, is known as a Cowell integration scheme. The scheme doesn't really care what kind of orbit it is generating, and it can handle deviations from classical two-body motion. Moreover, the user has complete control concerning its accuracy with his choice of propagation step-size. Figure 5-1 illustrates the common logarithm of the relative error in Shanks' 8-12 numerical integration scheme. The error is based on a 200 nm unperturbed orbit. In this instance, the Earth is modeled as a sphere so that, at least theoretically, the altitude is a constant 200 nm. The computed altitude will, of course, differ from the true altitude depending on the integration step-size and the number of integration steps. The figure illustrates the error as a function of the final time in days of 86400 seconds. (Equivalent to the duration of the orbit since the initial time was zero.) The relative error is shown because of its connection with the number of significant figures and independence from any particular units employed. From elementary numerical analysis, it is known that if the relative error in any number is not greater than 0.5 X 10-N, the number is certainly correct to N significant figures. It is seen, for example, that if an orbit is propagated with a step size of 300 seconds (5 minutes) for a duration of 7 days (2016 steps), the logarithm of the relative error is about -5.5. Since log(0.5 X 10-N) ≈ -(0.301 + N) and -5.5 < -5.301, the computation is accurate to 5 significant figures, which is equivalent to about 0.01 nm for a 200.00 nm orbit. It is emphasized that the results obtained from the figure are based on an unperturbed orbit. However, they will generally be valid for gravity perturbed orbits as well. Trajectory Maker Integration Scheme Algorithm Supplement 5-5 ────────────────────────────────────────────────────────── Things are quite different for objects entering the Earth's atmosphere. Such an object experiences a violent transient in its acceleration profile that can tax most any integrator's numerical integrity. Several factors are involved. The speed of the object is squared in the drag force equation, the density varies by orders of magnitude becoming large in the 100 km region and the drag force is sensitive to the ballistic coeffi- cient. (Small ballistic coefficients result in larger drag forces than large ballistic coefficients.) Trajectory Maker Integration Scheme Algorithm Supplement 5-6 ────────────────────────────────────────────────────────── Figure 5-1. Numerical Integration Error Trajectory Maker Integration Scheme Algorithm Supplement 5-7 ────────────────────────────────────────────────────────── References 1. Shanks, B. F., Solutions of Differential Equations by Evaluations of Functions, Math. Comp. 20 (1966) pp. 21-38. Trajectory Maker Coordinate Systems Algorithm Supplement 6-1 ────────────────────────────────────────────────────────── 6. COORDINATE SYSTEMS AND TRANSFORMATIONS We have seen that in order to propagate the equations of motion with a Cowell integration scheme, all that is required is an initial state vector given in Earth- centered-inertial coordinates, . . . [x, y, z, x, y, z]T at some epoch t0. Rectangular coordinates such as these are difficult to visualize. Nevertheless, the program should include the option for specifying input in the form of ECI coordinates, because these coordinates will be known for many applications; in particular for ballistic missile trajectories. However, there exist other coordinate systems which will simplify our visualization process when it comes to describing an initial trajectory with input data, and interpreting various characteristics with output data. 6.1 Time Considerations A fundamental parameter that is needed to define an inertial reference system is the Greenwich hour angle of the vernal equinox at epoch or orbit insertion. This angle is equivalent to Greenwich sidereal time at epoch and serves to fix a satellite ground trace relative to the rotating earth. We suppose that the initial time, or epoch, is specified by the year, I, month, J, day, K, and universal time in hours, minutes, and seconds (past midnight). These units may be converted to the Julian date JD, measured from January 0.0 by the clever expression in Wertz [4], JD = K - 32075 + 1461*( I + 4800 + (J - 14)/12)/4 + 367*( J - 2 - (J - 14)/12*12)/12 - 3*((I + 4900 + (J - 14)/12)/100)/4 - 0.50 where I, J, K are defined to be integers and the arithmetic is performed in integer mode. That is, floating point numbers are truncated at each stage of a computation, such as done by a computer language such as C or FORTRAN. Trajectory Maker Coordinate Systems Algorithm Supplement 6-2 ────────────────────────────────────────────────────────── The right ascension of the Sun is obtained using Newcomb's equation [1] Ru = 279.6910 + 36000.7689 T + 0.0004 T2 (degrees) where T is the number of Julian centuries of 36525 days elapsed since noon Greenwich mean time on January 0, 1900. This value is computed from T = (JD - JD0)/36525 where JD0 is the Julian day from Greenwich mean noon on January 0.5, 1900 which is 2,415,020.0. Universal time in degrees is UT = 15*(hours + minutes/60 + seconds/3600) So that the Greenwich sidereal time is GST = Ru + UT - 180 (degrees) 6.2 Orbital Elements A convenient method for defining an orbit is with its orbital elements. These elements describe its size, shape and orientation. They are: semimajor axis a, eccentri- city e, inclination i, longitude of the ascending node Ω, argument of perigee w, and some time when the orbiting object passes perigee τ. This latter element is also equivalent to mean anomaly M, but when specifying non-elliptical orbits, it is more convenient to specify time past perigee and compute mean anomaly from the formula M = nτ where n is the mean motion given by n = (µ a3)1/2 This parameter is the same as 2π/P for elliptic orbits having period P. It is a known from two-body theory that orbital elements at an instant in time, in particular at the epoch, are Trajectory Maker Coordinate Systems Algorithm Supplement 6-3 ────────────────────────────────────────────────────────── integration constants of the differential equations of motion, and these constants are equivalent to the instantaneous rectangular coordinates by a given coordinate transformation. This transformation is from Roy [3]. A geocentric, inertial coordinate system is defined in the orbit plane with principle axis x* along the major axis positive in the direction of perigee, another axis z* normal to the orbit plane in the direction of positive angular momentum and a third axis y* defined so that the system is right handed. From the orientation angles Ω, w, i, the direction cosines of the x*, y*, z* two planar axes relative to the ECI system may be obtained. They are l1 = cos Ω cos w - sin Ω sin w cos i m1 = sin Ω cos w + cos Ω sin w cos i n1 = sin w sin i l2 = -cos Ω cos w - sin Ω sin w cos i m2 = -sin Ω cos w + cos Ω sin w cos i n2 = cos w sin i These cosines are independent of the specific type of orbit be it elliptic, hyperbolic or parabolic. Elliptical Case ─────────────── The ECI state vector for the elliptical case is given by x = al1 cos E + bl2 sin E - ael1 y = am1 cos E + bm2 sin E - aem1 z = an1 cos E + bn2 sin E - aen1 . x = (bl2 cos E - al1 sin E)(na/r) . y = (bm2 cos E - am1 sin E)(na/r) . z = (bn2 cos E - an1 sin E)(na/r) where b is the orbit's semiminor axis given by Trajectory Maker Coordinate Systems Algorithm Supplement 6-4 ────────────────────────────────────────────────────────── b = a(1 - e2)1/2 and r is objects' radius vector magnitude given by r = a(1 - e cos E) The parameter E is the object's eccentric anomaly obtainable from Kepler equation E - e sin E = M by employing Newton's method for finding roots. Hyperbolic Case ─────────────── The ECI state vector for the hyperbolic case is given by x = -al1 cosh F + bl2 sinh F + ael1 y = -am1 cosh F + bm2 sinh F + aem1 z = -an1 cosh F + bn2 sinh F + aen1 . x = (bl2 cosh F - al1 sinh F)(na/r) . y = (bm2 cosh F - am1 sinh F)(na/r) . z = (bn2 cosh F - an1 sinh F)(na/r) where b is the orbit's semiminor axis given by b = a(e2 - 1)1/2 and r is objects' radius vector magnitude given by r = a(F cosh F - 1) The parameter F is the objects' hyperbolic areal parameter obtainable from the hyperbolic analog of Kepler's equation F sinh F - F = M also, by employing Newton's method for finding roots. Trajectory Maker Coordinate Systems Algorithm Supplement 6-5 ────────────────────────────────────────────────────────── Parabolic Case ─────────────── The ECI state vector for the parabolic case is given by x = l1r x* + l2r y* y = m1r x* + m2r y* z = n1r x* + n2r y* . x = l1 x* + l2 y* . y = m1 x* + m2 y* . z = n1 x* + n2 y* where x* = √(µp) ( sin f* cos f* /p - sin f* /r) y* = √(µp) ( sin f* sin f* /p - cos f* /r) For the parabolic case, the semimajor axis a is defined, herein, as the object's perigee distance and the conic parameter is computed from p = 2a The parameter f* is the object's true anomaly and may be found by solving Barker's cubic equation _ D + 1/3 D3 = 2n(t - τ) where D = tan (f*/2) and _ n2 p3 = µ Details of the solution are given in Roy [3]. Trajectory Maker Coordinate Systems Algorithm Supplement 6-6 ────────────────────────────────────────────────────────── 6.3 Geospherical Inertial Coordinates Geospherical inertial coordinates are spherical coordinates having their origin at the Earth's center. These coordinates consist of an objects range r, right ascension α and declination δ. Declination is equivalent to geocentric latitude, until now, denoted by φ'. These coordinates are defined by r = √(x2 + y2 + z2) δ = arcsin(z/r) α = arctan(y/x) In addition, the following flight parameters are computed: Inertial velocity . . . v = √(x2 + y2 + z2) Flight path angle Γ = arcsin(x'/v) Velocity azimuth angle σ = arctan(y'/z') where ┌─ ─┐ ┌─ ─┐ ┌─ ─┐ │ x' │ │ cos δ cos α cos δ sin α sin δ │ │ x │ │ │ │ │ │ │ │ y' │ = │ -sin α cos α 0 │ │ y │ │ │ │ │ │ │ │ z' │ │ -sin δ cos α -sin δ sin α cos δ │ │ z │ └─ ─┘ └─ ─┘ └─ ─┘ All the remaining coordinates used for the software system are relative to the rotating Earth. 6.4 Earth-Fixed-Geocentric Coordinates Earth-fixed-geocentric (EFG) coordinates have origin at Trajectory Maker Coordinate Systems Algorithm Supplement 6-7 ────────────────────────────────────────────────────────── the Earth's geocenter, fundamental plane coincident with the equatorial plane and principle axis in the direction of the Greenwich meridian. This coordinate system, being fixed relative to the Earth is actually a rotating coordinate system. Letting e, f, g denote these coordinates, the e-axis lies in the equatorial plane and is positive in the direction of Greenwich; the g-axis is coincident with the Earth's polar axis, i.e., its spin axis; and the f-axis lying in the equatorial plane such that the system is right-handed. This coordinate system is related to the inertial frame by the transformation e = x cos(gst) + y sin(gst) f = -x sin(gst) + y cos(gst) g = z where gst is the current Greenwich sidereal time. If gha is the Greenwich sidereal time at epoch, then gst= gha + Ωt (mod 360) where t is the time past epoch. 6.5 Local Topocentric Coordinates Local topocentric coordinates (UVW) have origin at some local site on the surface of the Earth, and fundamental plane tangent to the surface at the site. This reference system is also a rotating coordinate system. Letting u, v, w denote these coordinates, where the w-axis is normal to the tangent plane at the site, positive upwards; the principle u-axis is lies in the tangent plane, positive at in an azimuthal direction ; and the v-axis lies in the tangent plane such that the system is right-handed. It is shown in O'Connor [2], that the transformation from earth-fixed-geocentric coordinates to local topocentric coordinates is given by the transformation Trajectory Maker Coordinate Systems Algorithm Supplement 6-8 ────────────────────────────────────────────────────────── ┌─ ─┐ ┌─ ─┐ │ u │ │ e - e0 │ │ v │ = ABC│ f - f0 │ │ w │ │ g - g0 │ └─ ─┘ └─ ─┘ where the rotation matrices are ┌─ ─┐ │ sin cos 0 │ A = │ -cos sin 0 │ │ 0 0 1 │ └─ ─┘ ┌─ ─┐ │ 1 0 0 │ B = │ 0 sin φ cos φ │ │ 0 -cos φ sin φ │ └─ ─┘ ┌─ ─┐ │ -sin cos 0 │ C = │ -cos -sin 0 │ │ 0 0 1 │ └─ ─┘ where φ, denote the latitude and longitude of the site, respectively, and e0, f0, g0, are its EFG coordinates. These latter coordinates may be computed directly from ┌─ ─┐ │ a │ e0 = │ ──────────────── + h │ cos φ cos * │ √(1 - e2 sin2 φ) │ └─ ─┘ ┌─ ─┐ │ a │ f0 = │ ──────────────── + h │ cos φ sin * │ √(1 - e2 sin2 φ) │ └─ ─┘ ┌─ ─┐ │ a(1 - e2) │ g0 = │ ──────────────── + h │ sin φ │ √(1 - e2 sin2 φ) │ └─ ─┘ Trajectory Maker Coordinate Systems Algorithm Supplement 6-9 ────────────────────────────────────────────────────────── If the azimuth angle is 90°, this system is referred to as the East-North-Up coordinate system. In this case the azimuthal matrix is the unit matrix ┌─ ─┐ │ 1 0 0 │ A = │ 0 1 0 │ │ 0 0 1 │ └─ ─┘ and the transformation becomes ┌─ ─┐ ┌─ ─┐ │ ue │ │ e - e0 │ │ vn │ = BC│ f - f0 │ │ wu │ │ g - g0 │ └─ ─┘ └─ ─┘ where ┌─ ─┐ │ -sin cos 0 │ BC= │ -sin φ cos -sin φ sin 0 │ │ cos φ cos cos φ sin sin φ │ └─ ─┘ By differentiating and solving for the e, f, g velocity components, we find that . . . . . . [ e f g ]T = (BC)T [ u v w ]T These components may be identified with the wind components in Section 4. That is R= (BC)T so that ┌ ─┐ ┌─ ─┐ │ q1 │ │ qe │ │ q2 │ = R │ qn │ │ q3 │ │ qu │ └ ─┘ └─ ─┘ Carrying out the matrix operations, we have q1 = -[qe sin + qn cos * sin φ] q2 = [qe cos - qn sin * sin φ] q3 = qn cos φ Trajectory Maker Coordinate Systems Algorithm Supplement 6-10 ────────────────────────────────────────────────────────── which are the inertial wind velocity components required for the equations of motion. 6.6 Common Radar Coordinates Radar coordinates (RAE) are used for measurements of at some particular geographic location or site. They consist of an object's range, azimuth and elevation. Range R is the distance from the site to the object; azimuth A is the angle between the object and the radar site's meridian, measured clockwise from the northern direction; and elevation E is the angle above or below the horizon. These coordinates are given by the following expressions R = √(ue2 +vn2 + wu2) A = arctan(ue/vn) E = arcsin(wn/R) where the subscripts e, n, u denote the east, north and up components, respectively. In addition range rate is obtained by differentiating range. We have, . . . . R = (ueue + vnvn + wuwu)/R 6.7 Geodetic Coordinates and Computation of Altitude Altitude of an orbiting object is its height above the Earth's surface. If the Earth were a sphere of radius of radius a, then the altitude of an object at a distance r from the geocenter would be h = r - a But the equatorial and polar radii differ by nearly 12 nm; The polar radius being the smaller of the two. A parameter commonly used for measuring the shape of the Earth is its flattening defined by f = (a - b)/a where, in this section, a and b denote the Earth's Trajectory Maker Coordinate Systems Algorithm Supplement 6-11 ────────────────────────────────────────────────────────── equatorial and polar radii, respectively. Flattening is usually provided as a reciprocal, for example the value used herein, is the WGS-72 value of 1/298.26. An improved altitude model is h = r - a(1 - f sin φ') where, as before φ', denotes the declination of the object. The declination is equivalent to geocentric latitude. Both of these expressions are used by the trajectory software, but are principally used for estimates in the computer graphics. A great deal of effort was expended in order to accurately model the density profile, so we should expect an accurate altitude model from which density is obtained during orbit propagation. An altitude model that will achieve any desired degree of precision that is also extremely fast; indeed the one used in this paper, is an unpublished algorithm conceived by L. H. Ivy. A by-product of his algorithm are the trigono- metric functions required to construct the geocentric to topocentric rotation matrix, R sited in previous sections. Consider an instantaneous cross-section of the Earth defined by the plane containing its polar axis and the mass m. Define a two-dimensional coordinate system with u-axis formed by the intersection of this plane with the equatorial plane, and v-axis coincident with the polar axis. The cross-section is bounded by an ellipse whose equation is u2 v2 ── + ── = 1 a2 b2 Since the eccentricity of the ellipse is given by b2 = a2(1 - e2) the its equation may be written Trajectory Maker Coordinate Systems Algorithm Supplement 6-12 ────────────────────────────────────────────────────────── u2(1 - e2) + v2 = a2(1 - e2) Now, the tangent of the geodetic latitude φ of a point on the ellipse is the negative reciprocal of the slope of the line tangent to the ellipse at that point. It follows by differentiation that 1 v tan φ = ──────── ─── (1 - e2) u Consequently, a cos φ u = ──────────────── √(1 - e2 sin2 φ) a(1 - e2) sin φ v = ──────────────── √(1 - e2 sin2 φ) These coordinates are the coordinates of a point on the ellipse whose geodetic latitude is φ. On the other hand, the coordinates of a point at an altitude h whose geodetic latitude is φ are uh= u + h cos φ vh= v + h sin φ The geocentric latitude φ' of this point is given by ┌─ ─┐ │ e2 │ tan φ' = │ 1 - ─────────────────────────── │ tan φ │ 1 + (h/a) √(1 - e2 sin2 φ ) │ └─ ─┘ One can also find the following expression for altitude h= uh cos φ + vh sin φ - a√(1 - e2 sin2 φ) These two equations form the basis of Ivy's algorithm. In essence, the algorithm is an iterative process which uses these equations starting with geocentric latitude as an estimate for geodetic latitude. Trajectory Maker Coordinate Systems Algorithm Supplement 6-13 ────────────────────────────────────────────────────────── To complete the algorithm, Ivy defines a miss distance as the distance between the current estimate for geodetic latitude φi and the line normal to the reference ellipsoid corresponding to the true geodetic coordinates. Let ui, vi be a point on the ellipsoid whose latitude is φi. Then the equation of the line normal to the ellipsoid at this point, expressed in normal form, is (v - vi)cos φi - (u - ui)sin φi = 0 The distance from the line to the point u, v is p = │ (u - ui)sin φi - (v - vi)cos φi │ Using (3) and (4) with u = ui, v = vi, we have the convergence criterion, │ ae2sin φi cos φi │ p = │ uh sin φi - vh cos φi - ───────────────── │ │ √(1 - e2 sin2 φ) │ The algorithm can now be summarized. Algorithm Summary ───────────────── 1. Input the coordinates of the object: x, y, z 2. Compute: √(x2 + y2) 3. Compute: tan φ' = z/√(x2 + y2) 4. Replace: tan φ = tan φ' 5. Compute: cos2 φ = 1/(1 + tan2 φ) 6. Compute: sin2 φ = 1 - cos2 φ 7. Compute: rdcl = √(1 - e2 sin2 φ) 8. Compute: cos φ = √cos2 φ 9. Compute: sin φ = cos φ tan φ Trajectory Maker Coordinate Systems Algorithm Supplement 6-14 ────────────────────────────────────────────────────────── 10. Compute: h = x cos φ + y sin φ - a rdcl 11. Compute: p = │x sin φ - y cos φ - ae2sin φ cos φ/rdcl│ 12. Test: If p is less than some error criterion, end algorithm. Otherwise, 13. Compute: D = 1 + (h/a) rdcl 14. Compute: tan φ = tan φ '/(1 - e2/D) 15. Repeat Steps 5 through Step 14. Optional computations: Compute geocentric and geodetic latitudes: φ' = arctan(tan φ'), φ = arctan(tan φ) Compute right ascension: cos α = x/√(x2 + y2), sin α = y/√(x2 + y2) α = arctan(sin α /cos α) It is worthy to note that no trigonometric functions are actually computed with this algorithm, yet it makes available the sine and cosine of the geodetic latitude. Thus the rotation matrix, required for computing wind velocity components, may readily be constructed. This is a great savings in computer resources because this matrix must be computed for each integration step. The algorithm was tested with latitudes ranging from -90° to 90°, and altitude ranging from -300 to 5 X 1012 feet at each latitude. Miss distances of 10-1, 10-4 and 10-5 feet were used as convergence criteria. In no case were more than three iterations required, and no more than two were required for miss distances of 10-4 feet. Trajectory Maker Coordinate Systems Algorithm Supplement 6-15 ────────────────────────────────────────────────────────── References 1. Anon., United States Navel Observatory and Royal Greenwich Observatory, Explanatory Supplement to the American Ephemeris and Nautical Almanac, Her Majesty's Stationary Office, London, 1961. 2. O'Conner, J. J., Transformations Applicable to Missile and Satellite Trajectory Computations, RCA International Service Corporation, Technical Report AFETR-TR-75_29, Patrick Air Force Base, Florida, July 1975. 3. Roy,A. E., The Foundations of Astrodynamics, The Macmillan Company, New York, 1965. 4. Wertz, J. R., ed., Spacecraft Attitude Determination and Control, D. Reidel Publishing Company, Boston, 1985. Trajectory Maker Bibliography Algorithms Supplement 7-1 ────────────────────────────────────────────────────────── 7. BIBLIOGRAPHY Abell, G., Exploration of the Universe, Holt, Reinhart and Winston, New York, 1964. Anon., Space Flight Handbooks, Volume I, NASA SP-33 Part 1, Office of Scientific and Technical Information, NASA, Washington D. C., 1963. Anon., United States Navel Observatory and Royal Greenwich observatory, Explanatory Supplement to the American Ephemeris and Nautical Almanac, Her Majesty's Stationary Office, London, 1961. Anon., Orbital Flight Handbook, Office of Scientific and Technical Information, National Aeronautics and Space Administration, Washington D.C., 1963. Bates, R., and D. Mueller, J. White, Fundamentals of Astrodynamics, Dover Publications, Inc., New York, 1971. Battin, B., An Introduction to the Mathematics and Methods of Astrodynamics, American Institute of Aeronautics and Astronautics, Inc., New York, 1987. Deetz, C. and O. Adams, Elements of Map Projections, Greenwood Press, New York, 1969. Felberg, F., New One-Step Integration Methods of High- -Order Accuracy Applied to Some Problems in Celestial Mechanics, NASA Technical Report NASA TR R-248, George C. Marshall Space Flight Center, Huntsville, Alabama, October 1966, pp. 32-33. Gaposchkin, F. M., ed., 1973 Smithsonian Standard Earth (III), Smithsonian Astrophysical Observatory, Special Report 353, Cambridge, Massachusetts, November 28, 1973. Hieb, H., Western Test Range Geodetic Coordinate Manual, Part 1, Federal Electric Corporation, Feb. 1980. Ivy L. and H. Reynolds, On the Sensitivity of Impact Prediction to Methods of Interpolating Atmospheric Density, ITT Technical Note, TN-76-1463, Vandenberg Air Force Base, 1976. Maling, D., Coordinate Systems and Map Projections, George Phillip and Son Limited, London, 1973. Trajectory Maker Bibliography Algorithms Supplement 7-2 ────────────────────────────────────────────────────────── Mc Cuskey, S., Introduction to Celestial Mechanics, Addison Wesley Publishing Company, Inc., Palo Alto, 1963. Menzner, et. al., Defining Constants, Equations, and Abbreviated Tables of the 1975 U. S. Standard Atmosphere, National Aeronautics and Space Administration, Technical Report NASA TR R-459, Washington D. C., May 1976. O'Connor, J. J., Transformations Applicable to Missile and Satellite Trajectory Computations, RCA International Service Corporation, Technical Report AFETR-TR-75-29, Patrick Air Force Base, Florida, July 1975. Pearson, F., Map Projections: Theory and Applications, CRC Press, Inc., Boca Raton, 1990. Ratcliffe, J., A., An Introduction to the Ionosphere and Magnetosphere, Cambridge Univ. Press., 1972, p.31. Ridpath, I., Longman Illustrated Dictionary of Astronomy and Astronautics, York Press, Essex, 1988. Roy, R., The Foundations of Astrodynamics, The Macmillan Company, New York, 1965. Shanks, B. F., Solutions of Differential Equations by evaluations of Functions, Math. Comp. 20 (1966) pp. 21-38. Steers, D., An Introduction to the Study of Map Projections, University of London Press Ltd., London, 1962. Struick, D., Lectures on Classical Differential Geometry, Addison-Wesley Publishing Company, Inc., Reading, 1961. Thompson, R. F., Spline Interpolation on a Digital Computer, Goddard Space Flight Center, Technical Report X-692-70-261, Greenbelt, Md., July 1970. Waison, G. R., Prediction of Orbital Decay Rates Due to Atmospheric Drag, TRW Program Technical Report 3150-6020- RO-000, January 7, 1966. Wertz, J. R., ed., Spacecraft Attitude Determination and Control, D. Reidel Publishing Company, Boston, 1985. Whittaker, F. and G. Watson, A Course in Modern Analysis, Cambridge University Press, Fourth Ed. 1963.